Problem 1. (60 points)

Let us fix a natural number q. Set n = q^2. The matrix A will be the n \times n matrix with entries¹

A_{ij} = \begin{cases} \gamma, & i=j,\\ -1, & \Big(|i - j| = 1 \text{ and } \min(i,j) \mod q \neq 0\Big) \text{ or } |i - j| = q,\\ 0, & \text{otherwise}. \end{cases}

Recall that \min(i,j) \mod q \neq 0 means that \min(i,j) is not divisible by q.

This is how the matrix looks for q = 2 and \gamma = 4:

\begin{pmatrix} 4&-1&-1& 0\\ -1& 4& 0&-1\\ -1& 0& 4&-1\\ 0&-1&-1& 4 \end{pmatrix} and q = 3: \begin{pmatrix} 4&-1& 0&-1& 0& 0& 0& 0& 0\\ -1& 4&-1& 0&-1& 0& 0& 0& 0\\ 0&-1& 4& 0& 0&-1& 0& 0& 0\\ -1& 0& 0& 4&-1& 0&-1& 0& 0\\ 0&-1& 0&-1& 4&-1& 0&-1& 0\\ 0& 0&-1& 0&-1& 4& 0& 0&-1\\ 0& 0& 0&-1& 0& 0& 4&-1& 0\\ 0& 0& 0& 0&-1& 0&-1& 4&-1\\ 0& 0& 0& 0& 0&-1& 0&-1& 4 \end{pmatrix}

It is probably easier to see the distribution of the values in the matrix plot for q = 4:

Note that every qth element on the diagonals next to the main diagonal is 0.

Do the following. Make sure that the style of your plots is similar to the example solutions below!

1. Modify the preconditioned CG method program from Lecture b3 with the above matrix. The program should read q, \gamma and \omega from the keyboard input (in this order). Use {\varepsilon}= 10^{-10} and b_i = 1/ (q +1)^2. The program should print the norm of the residual \|r\| for after every iteration.

   Submit the code.

2. (20 points) Run the program for q = 50, \gamma = 4 and each \omega = \{1., 1.1, 1.2, \ldots, 1.9\}, find the number of iterations that the program requires for given \omega. Make a plot of the number of iterations as a function of \omega in gnuplot. Use with linespoints.

   Submit the plot with your student ID as the title.

3. (20 points) Run the program for q = 50, \gamma = 5 and each \omega = \{1., 1.1, 1.2, \ldots, 1.9\}, find the number of iterations that the program requires for given \omega. Make a plot of the number of iterations as a function of \omega in gnuplot. Use with linespoints.

   Submit the plot with your student ID as the title.

4. Modify the program so that it does not use any preconditioning (you can just have function precond return z = r). Find the number of iterations the program requires for q = 50, \gamma = 4.

   Submit the number of iterations.

Problem 2. (40 points)

Modify the program from Problem 1 to print the values of the solution x in the following format:

Example: This is the output for q = 4 and \gamma = 4:

3.33333313E-02   4.66666631E-02   4.66666631E-02   3.33333313E-02
4.66666631E-02   6.66666627E-02   6.66666627E-02   4.66666631E-02
4.66666631E-02   6.66666627E-02   6.66666627E-02   4.66666631E-02
3.33333313E-02   4.66666631E-02   4.66666631E-02   3.33333313E-02

You can use the following code to print the values in the above format:

do i = 1, q
    print *, x((i-1) * q + 1:i * q)
enddo

If you store the output in a file (say data.txt), you can make a 3D plot with contours in gnuplot using the following commands

set title '123456789'
unset key
set contour base
set cntrparam levels 50
splot 'data.txt' matrix with lines

Always use {\varepsilon}= 10^{-10}.

Do the following. Make sure that the style of your plots is similar to the example solutions below!

1. Run the program with q = 50, \gamma = 4 and \omega = 1.8, and b_i = 1/(q + 1)^2. Plot the solution x as explained above.

   Submit the plot with your student ID as the title.

2. Run the program with q = 50, \gamma = 4 and \omega = 1.8, and b_i = \begin{cases} 1, & i = q(q + 1) / 2,\\ 0, & \text{otherwise}. \end{cases} Plot the solution x as explained above.

   Submit the plot with your student ID as the title.

3. Run the program with q = 60, \gamma = 4 and \omega = 1.8, and b_i = \begin{cases} 1, & i = q(q+1)/3,\\ -1, & i = 2 q(q+1) / 3,\\ 0, & \text{otherwise}. \end{cases} Plot the solution x as explained above.

   Submit the plot with your student ID as the title.

4. Run the program with q = 60, \gamma = 4 and \omega = 1.8, and b_i = \begin{cases} \sin(4\pi i /(q+1)), & i \leq q,\\ 0, & \text{otherwise}. \end{cases} Plot the solution x as explained above.

   Submit the plot with your student ID as the title.